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Abstract 

Starting from the adiabatic time-dependent Hartree-Fock approximation (ATDHF), we propose 
an efficient method to calculate the Thouless-Valatin moments of inertia for the nuclear system. 
The method is based on the rapid convergence of the expansion of the inertia matrix. The accuracy 
of the proposed method is verified in the rotational case by comparing the results with the exact 
Thouless-Valatin moments of inertia calculated using the self-consistent cranking model. The 
proposed method is computationally much more efficient than the full ATDHF calculation, yet it 
retains a high accuracy of the order of 1%. 

PACS numbers: 21.60.Jz, 21.60.Ev, 21.10.Re 
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I. INTRODUCTION 



The variation of nuclear ground-state shapes is governed by the modification of the shell- 
structure of single-nucleon orbitals. Far from the valley of ^-stability, the energy spacings 
between single-nucleon levels change considerably with the number of neutrons and/or pro- 
tons. The reduction of spherical shell closure is often associated with the occurrence of 
deformed ground states and, in many cases, with the phenomenon of coexistence of different 
shapes in a single nucleus. A quantitative description of the evolution of nuclear shapes, 
including regions of short-lived exotic nuclei that are becoming accessible in experiments 
at radioactive-beam facilities, necessitate accurate modeling of the underlying microscopic 
nucleonic dynamics. Major advances in nuclear theory have recently been made in stud- 
ies of complex shapes and the corresponding excitation spectra and electromagnetic decay 
patterns, especially in the framework of nuclear energy density functionals (EDFs) jl-5]. 

A microscopic, EDF-based description of complex collective excitation spectra usually 
starts from a constrained Hartree-Fock plus BCS (HFBCS) or Hartree-Fock-Bogoliubov 
(HFB) calculation of the binding energy surface with the mass multipole moments as con- 
strained quantities. The static nuclear mean-field is characterized by symmetry breaking: 
translational, rotational and particle number. Even though symmetry breaking incorporates 
important static correlations (e.g., deformations and pairing), the static self-consistent so- 
lution can only provide an approximate description of bulk ground-state properties such as 
masses and radii. Modeling excitation spectra and transition rates in the EDF framework 
necessitates a systematic treatment of dynamical effects related to restoration of broken 
symmetries and fluctuations in collective coordinates. 

One possible approach to five- dimensional quadrupole dynamics that restores rotational 
symmetry and allows for fluctuations around triaxial mean-field minima is to formulate 
a collective Hamiltonian, with deformation-dependent inertia parameters determined by 
microscopic self-consistent mean-field calculations. The dynamics of the collective Bohr 
Hamiltonian is governed by the vibrational inertial functions and the moments of inertia js]. 
For these quantities either the Gaussian overlap approximation of the generator coordi- 
nate method (GCM-GOA) (Yoccoz masses [3]) or the adiabatic time-dependent Hartree- 
Fock-Bogoliubov (ATDHFB) expressions (Thouless-Valatin masses [8|) can be used. The 
Thouless-Valatin masses have the advantage that they also include the time-odd components 
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of the self-consistent mean field and, in this sense, the full dynamics of a nuclear system. 
This can be seen most clearly in the case of translational motion, where the Thouless-Valatin 



mass corresponds to the exact mass A ■ m of A nucleons jo], whereas the GCM method pro- 
duces the exact value only when the center of the mass velocity is also included as the 



generator coordinate [10[ . The calculation of the Thouless-Valatin masses is often simplified 
by adopting the cranking formulas [lj], [3] that neglect the residual interaction. In that 
case the Thouless-Valatin corrections are usually taken into account by scaling the inertia 



parameters with an empirical factor (w 1.2 - 1.4) |l3l-ll5l]. 

In this work we present an efficient method to calculate the Thouless-Valatin moments of 
inertia for the nuclear system. The method is based on the rapid convergence of the expan- 
sion of the inertia matrix. The accuracy of the proposed method is verified in the rotational 
case by comparing the results with the exact Thouless-Valatin moments of inertia calculated 
using the self-consistent cranking model. The proposed method is computationally much 
more efficient than the full ATDHF calculation, yet it retains a high accuracy of the order 
of 1%. 

II. THEORETICAL FRAMEWORK 

We begin with a brief review of the adiabatic time-dependent Hartree-Fock theory. A 



more detailed exposition of this formalism can be found, for instance, in Refs. [ji], 



3- The 



aim of the ATDHF theory is to derive in a fully microscopic and consistent way a Hamiltonian 
for the description of collective phenomena in which many nucleons act coherently. The 
theory is based on two approximations: i) in the TDHF one assumes that the many-body 
time-dependent wave function of the system is a Slater determinant at all times; and ii) in 
the adiabatic approximation the collective motion is slow compared to single-particle motion 
and, therefore, the collective kinetic energy is a quadratic function of the velocities. 

To identify the components of the density matrix that correspond to the coordinates 
and momenta of the collective Hamiltonian, we recall that the coordinates are even and the 
momenta are odd under time-reversal, and decompose the density matrix in the following 
way: 

p(t) = e ix{t)/h p (t)e- ix{t)/h . (1) 
Both matrices, po(t) and are Hermitian and time-even. p (t) represents the coordinates 
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of the collective Hamiltonian, and x(t) is the "adiabaticity parameter" that must be small 
compared to unity. At all times p(t) is a Slater determinant, that is, po{t) 2 = po(t) and 
Trpo = N, N being the particle number. In the following we work in the basis in which p is 
diagonal, and consequently the operators po an d <j = 1 — p project onto hole and particle 
states, respectively. This basis depends on time because po(t) is a function of time. 

In the adiabatic approximation it is assumed that the total density p(t) of the system 
is always close to the density po(t), that is, the matrix \ that introduces the time-odd 
components remains small at all times. Expanding the density matrix to second order in 
the operator x, the following expression is obtained: 

p "{ l + \ x ~ h x ") p ° ( x - 5* - h x ") 53 + Pi + fe (2) 

where 

i i 

Pi = K [X, Po] = ^ (xpo ~ PoX) , (3) 

1 11 

P2 = [x, [X, Po\] = j^XPoX - (x 2 Po + PoX 2 ) ■ (4) 

Pi is linear in x, time-odd, and has only ph and hp non-vanishing matrix elements. p2 is 
quadratic in x, therefore time-even, and has only hh and pp matrix elements. The many- 
body Hamiltonian can also be expanded to second order in the operator x- 

hab(p) = tab + V a dbcPcd = tab + V a dbc(po)cd + V a dbc(Pl)cd + V a dbc{p2) cd, (5) 

cd cd cd cd 

where t is the kinetic energy operator, and V denotes a generic two-body interaction. The 
Hamiltonian contains time-even (/i and T 2 ) and time-odd parts (Ti) 

h(p) = h + T 1 + T 2 . (6) 

Consequently the time-dependent Hartree-Fock equation ihp = [h, p] also decomposes into 
two equations: 

ihpo = [ho, px] + \T U po], (7) 
ihp! = [h , po] + [Ti, pi] + [r 2 , po]- (8) 

In Eq. (jS|) the term [ho, pi] has been neglected because the ph and hp components are small, 



and the pp and hh parts vanish 



171 ] . The total energy of the system 

E = tabPba + 7} ^2 PbaVadbcPcd, (9) 
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ab abed 



can be expressed in terms of the variables po and p±, or po and x- Terms which depend on 
the velocity in second order build the kinetic energy of the collective Hamiltonian: 

K = ^ ( h o)ab(p2)ba + - ^2 (Pl)baV a dbc(pi)cd- (10) 
ab abed 

We recall that the matrix po projects onto the hole states and, inserting Eqs. ([3]) and (J4]) in 
the expression above, the kinetic energy can be written: 

(11) 

j \r j 

where the matrix A is Hermitian, and B is symmetric 

Aphp>h> = (e P — £h)5pp'$hh> + Vphh'p', Bphp'h' = Vphp'h'- (12) 




In Ref. 18| it is shown that the effective p/i-interaction in relativistic point coupling 



models can be written as a sum of separable terms 

V^ = £<&W3 (is) 

r 

The same is true for relativistic Hartree models with meson exchange forces. The single 
particle operators Q r are either even or odd under time- reversal. The time-odd operators 
correspond to the isoscalar and isovector currents j(r) and fj(r). Implementing the inter- 
action ffl3l) into Eq. (TTUj) we find 



~ £ (pi) b aV ad bc(piU = ^Tr(g> 1 )KTr(g»*. (14) 

abed r 

Since p\ is time-odd the traces vanish for time-even operators and only the time-odd oper- 
ators in the matrices A and B contribute to the inertia parameters. 
The equation of motion ([7j) can be written into the following form: 

pA i (a -b\ I , \ ,. / 




Po J \ -o A J \ X ) \X 

To perform realistic calculations the dimension of the problem has to be reduced, that 
is, one has to select a small number of active degrees of freedom qi,...,qf. This means 
that we are able to generate a subset of time-even Slater determinants, characterized by the 
parameters q, with the following property: the solution of the ATDHF problem will always 



remain within this subset of Slater determinants. In other words, we have found a path 
Po(t) = po[q(t)] from which we can calculate the velocity 

Po(t)=q(t)^- (16) 
Next we define the operator P with the relation: 

£ = (17) 

and obtain the following expression for the kinetic energy: 

/ 



1 f 

K = oYl M Mw» > ( 18 ) 



where M MM /(q) denotes the real collective mass tensor 



M w ,(q) = l(p.-P)^( 



(19) 



To evaluate M, we have to invert the matrix 



AT 1 = I A B I =7W~ 1 + V (20) 
-S* A* 

in Eq. ([15]) . For this purpose we decompose the matrix into a diagonal part containing the 
energies of particle and hole states 

(■^o 1 )^ = ( 6 p - e h) Spp'S w , (21) 

and the residual interaction V, and use the fact that the interaction matrix elements are in 
most cases much smaller than the p/i-energies. This is because only the time-odd compo- 
nents of the residual interaction contribute. Therefore the matrix M. can be written in the 
following form: 

M = [Mq 1 + = Moit + VMoV 1 (22) 
We expand the factor in the square bracket and obtain: 

M = M - MqVMq + MqVMqVMq H • (23) 
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Since M 1 is diagonal, inverting this matrix is trivial and the problem is reduced to simple 
matrix multiplications. The zero-order term, of course, yields the Inglis-Belyaev formula: 

AC(q)4(p' -p)m(_;| 4E^. (24) 

The first- and the second-order terms 

Mi = -M VM , M 2 = M VM VM . (25) 

represent the leading corrections to the Inglis-Belyaev formula. The purpose of this ex- 
ploratory study is to determine the convergence of the expansion (I23I) . as well as the level of 
agreement with the Thouless-Valatin formula. In this work we only consider the moments 
of inertia for collective rotation, that is, the operator P corresponds to the components of 
the angular momentum vector . 

In fact, for a stationary deformed solution without external constraint, as it is discussed 
in the following application, the RPA-equation has a Goldstone mode. As discussed in detail 
in Sect. 8.4.7 of Ref. , from rotational invariance of the Hamiltonian, i.e. from [H, J x ] = 
we obtain for P = J x a spurious solution 

P 





= (26) 
On the other side, from Eq. ( fl9l) we see that we have to solve the inhomogeneous equation 

(27) 





Such a solution exists, because the inhomogeneous part of this equation is orthogonal to 
the Goldstone mode of Eq. (126|) . Of course, the explicit inversion of the matrix M~ l is 
technically complicated because it has to be carried out in the space orthogonal to the 
Goldstone mode. However, the method proposed here avoids these technical complications. 
In each order of the approximation the matrix fl23|) acts on the vector 

(28) 

which eliminates all spurious contributions. We also have to emphasize, that the problem 
of the Goldstone mode occurs only at the stationary points of the energy surface, where 




the constraint vanishes. For all other solutions the constraining operator does not commute 
with the angular momentum and therefore there exist no spurious solution. 

As a specific example of the nuclear energy density functional we consider the point- 



coupling implementation of a relativistic EDF - the functional PC-F1 
Ermf = J dr£ RMF (r) = ^2 J drijj i (r){-i'yV + m)^?) 



i=i 

+ / dr 



+ ^ TV (j T v)^UTvT + ^tvUtv)^UtvT + -^OiTSPTS + ^TSpTS^pTS + \pV^ 

(29) 

where ip denotes the Dirac spinor field of a nucleon, and the local isoscalar and isovector 
densities and currents 

A 

ps(r) = 5^(rM(r), (30) 
i=i 

A 



p TS (r) = Y,Mr)r 3 ^(r), (31) 

8=1 

A 

j»{r)= J]^(r)7^(r), (32) 
i=i 

A 

AW = ^^W7^(r), (33) 



8=1 



are calculated in the no-sea approximation: the summation runs over all occupied states 
in the Fermi sea. This means that only occupied single-nucleon states with positive energy 
explicitly contribute to the nucleon self-energies. In Eq. ( 129]) p p is the proton density, and 
A denotes the Coulomb potential. 

The matrix elements of the residual interaction are derived from the EDF Eq. (12 9p 

d 2 E RMF 

Vadbc — q q ) 1<J 4 J 

OpbaOPcd 

where generic indices (a, b, c, d, . . . ) denote quantum numbers that specify the single-nucleon 
state {ifj a }- These belong to three distinct sets: the index p (particle) denotes unoccupied 
states above the Fermi sea, the index h (hole) is for occupied states in the Fermi sea, and 
with a we denote the unoccupied negative-energy states in the Dirac sea. The calculation 
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of the moments of inertia involves only the time-odd terms of the residual interaction, for 
which the isoscalar-vector channel plays the dominant role. The time-odd contributions of 
the isovector-vector and the electromagnetic fields are omitted because the corresponding 
couplings are small in comparison to the isoscalar-vector coupling. Here we make a further 
simplification by assuming that the nonlinear and the derivative terms can be neglected, 
that is, it is sufficient to retain only the linear isoscalar-vector term (see Fig. [TJand Tab. [I]): 

V adbc = -a v J [ipla^liplaip^r . (35) 
III. NUMERICAL TEST 
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FIG. 1: (Color online) The moment of inertia T x of 154 Sm computed using the cranking RMF 
framework. The various curves correspond to calculations that successively include the following 
terms of the residual interaction: the time-even part (Inglis-Belyaev formula denoted by IB), the 
linear (VL), nonlinear (VNL), and derivative (VD) time-odd terms in the isoscalar-vector channel. 



To verify our assumption that the time-odd part of the residual interaction can be ap- 
proximated by the linear vector term, we have analyzed the contributions of the different 
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TABLE I: Contributions to the moment of inertia X x of different terms of the residual interaction: 
the time-even part (Inglis-Belyaev formula denoted by IB), the linear (VL), nonlinear (VNL), 
and derivative (VD) time-odd terms in the isoscalar-vector channel. The calculation is performed 
using the cranking RMF framework with the PC-F1 interaction, and the cranking frequency is 
n = 0.001 MeV. 

T IB T VL -7-VNL T VD 
■^x x -^x -^x 

55.53 18.99 -0.31 1.44 



time-odd terms to the moments of inertia by performing a self-consistent cranking calcula- 
tion (see Refs. 3), 2o| and references cited therein). In the cranking framework there are 
two types of moments of inertia: the kinematic (or static) moment of inertia and the 

dynamic moment of inertia . They are defined as follows 

j^n) = i, j* (a) = %. (36) 

In a self-consistent calculation with very small vlaues of the rotational frequency, J^ 2 \VL) 
is identical to the Thouless-Valatin moment of inertia, the linear response to the external 
Coriolis field. At the band-head in even-even nuclei the two quantities and coincide 
and we use in the figures the character X for this quantity. Calculations that neglect the 
time-odd fields and take into account only the Coriolis operator QJ^ in the Dirac equation, 



underestimate the empirical moments of inertia by roughly 30% [21] . As an illustrative 
example, in Fig. [1] we plot the dynamic moment of inertia for the ground state band in 
154 Sm. By including only the linear time-odd term (VL) in isoscalar-vector channel the 
moment of inertia is enhanced by 34%, while the remaining two contributions: the non-linear 
term (VNL) and the derivative term (VD) yield less than 3%. The results are summarized 
in Tab. [H where we list the contributions of the linear vector, nonlinear vector and vector 
derivative terms to the moment of inertia. Thus in the remaining calculations the model 
includes only the linear time-odd term. 

To estimate the convergence of the expansion formula (12"B"1) , we have used it to calculate the 
moment of inertia I x for the 154 Sm isotope, in comparison with the exact Thouless-Valatin 
moment of inertia computed with the cranking code. In Fig. [2] the latter is compared with 
the zeroth-, first-, and second-order in the expansion Eq. (12"3~|) . Moreover, we also display the 
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FIG. 2: (Color online) The moment of inertia I x of 154 Sm computed at the zeroth-, first-, and 
second-order in the expansion Eq. (|23p based on the PC-F1 density functional, and compared with 
the values obtained at the corresponding iteration steps in the cranking RMF based on the same 
density functional. 

cranking results for each iteration step in the self-consistent cranking calculation, starting 
from the stationary solution without the cranking term (cf. Appendix). As one expects, the 
moment of inertia obtained after the first iteration is equal to the Inglis-Belyaev moment of 
inertia. The next two iterations are compared to the first and second order in the expansion 
formula (T23]) . We note that the values obtained after the second and third iteration steps are 
in complete agreement with the first- and the second-order corrections, respectively. In the 
Appendix we demonstrate that these values have to be identical, thus the results displayed 
in Fig. [2] provide a crucial test for the numerical implementation of the expansion Eq. (I23jl . 
Further iteration steps contribute to the value of the moment of inertia by less than 1%, 
that is, the convergence is quite rapid. We also emphasize that it is necessary to include the 
contributions from the negative-energy single-nucleon Dirac states in the calculation of the 
matrix Ai. Omitting the negative energy states leads to a significant overestimation of the 
second-order correction to the moment of inertia. 
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FIG. 3: (Color online) The moments of inertia I x of 152_164 Sm. The values computed at the zeroth- 
, first-, and second-order in the expansion Eq. (|23p are compared with the exact Thouless-Valatin 
moment of inertia calculated using the RMF cranking model (upper panel). For the same nuclei 
the ratio of the Thouless-Valatin and Inglis-Belyaev moments of inertia (lower panel). 

Finally, in Fig. [3] we display the moments of inertia X x for the sequence of even-even 
isotopes 152_164 Sm. The values computed at the zeroth-, first-, and second-order in the ex- 
pansion Eq. ( 1231) are compared with the exact Thouless-Valatin moment of inertia calculated 
using the RMF cranking model. Through the whole isotopic chain the expansion method, 
truncated to second order, yields values very close to the exact Thouless-Valatin moments 
of inertia, with the relative deviation ~ 1%. We note that the enhancement of the moment 
of inertia in comparison to the Inglis-Belyaev value ranges of 1.28 ~ 1.39. 
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IV. SUMMARY AND OUTLOOK 



Starting from the adiabatic time-dependent Hartree-Fock theory, we have introduced 
an efficient approximate method to calculate the Thouless-Valatin moments of inertia for 
the nuclear system. The method is based on the fact that the expansion of the inertial 
parameters converges rapidly because the matrix elements of the time-odd components of the 
residual interaction are usually small in comparison to the p/z-energies. This approximation 
is computationally much less demanding than the full ATDHF calculation, yet it retains high 
accuracy of the order of 1%. The accuracy of this method has been verified by comparing 
the results to the exact Thouless-Valatin rotational moments of inertia calculated within the 
cranking model. 

One might, of course, encounter problems in regions of level crossings, where the ph 
energies are no longer necessarily small compared to the matrix elements of the residual 
interaction V. In that case the matrix .M -1 in Eq. ( 120]) has to be decomposed in a different 
way as, for instance, by shifting the diagonal elements of V to Ai® 1 , or by adding and 
subtracting complex diagonal elements. 

The present study has been limited to the rotational moments of inertia. In future 
investigations we plan also the calculations of vibrational masses. In this case the momentum 



operator P is not known a priori. Several ways have been proposed in the literature 
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23[ to attack this problem. We hope to solve this problem by a similar expansion as 



in Eq. ( )23l) . Of course this method can also be used with different density functionals by 
simply replacing the time-odd residual interaction. Pairing correlations can be included 
by expanding the inverse of the QRPA matrix instead of the RPA matrix. Work in this 
direction is already in progress. 
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Appendix: Iterative solution of the cranking equation 

In this appendix it is demonstrated that the moment of inertia calculated at each step 
of the iterative solution of the cranking equation, coincides with the corresponding order of 
the expansion introduced in Sec. [Til We assume that the cranking frequency in the equation 
of motion 

[h(p)-nj x ,p] = 0, (37) 

is an infinitesimal quantity, that is, second and higher order terms in Q can be safely ne- 
glected. As the initial point we choose the self-consistent solution po for frequency Q = 0. 
The corresponding equation of motion reads: 

[ho, Po] = 0. (38) 

In the first step of the iteration we diagonalize the operator ho — £lj x , and compute the 
density po + Spo determined by the following relation: 

[h ,5po] = n[j Xl p ]. (39) 

In the basis which diagonalizes ho, the only no n- vanishing matrix elements of 5po are ph and 
hp. Using the definition of the matrix Ai Eq. (1211) . we obtain 

Sp °)=n M J } A. (40) 

In the following the shorthand notation is used: 

5p = I 6p ° ) and j x =\ jx ) , (41) 

that is, 5p~o = QA4oj x . After the first iteration we obtain the Inglis-Belyaev moment of 
inertia: 

h = = i ( a * ) ( Jj J = a*** s /,B < 42 > 
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In the second iteration we diagonalize the operator: 

hi - ttjx = h(p + 5p ) - Qj x = h + VSp - Qj x = h - Q. (j x - VM j x ) , (43) 
where VAiojx denotes the matrix 

(VM j x ) ab = ^ Vahb P M php'h'(jx)p>h> +VapbhMohph>p'{3x)* plhl ■ (44) 
php'h' 

The density p + 8 Pi is the solution of the equation of motion 

[h - n ( 3x - VM ] X ) , po + $Pi] = . (45) 

Again, 5p\ has only p/i-matrix elements and, therefore, we need only these elements of the 
matrix (]44p and find: 

5p 1 = QM (t-VMo)j x . (46) 
The moment of inertia obtained in the second iteration coincides with that defined by Eq. 

A25D 

h = ^Tr(jJ Pl ) = ( j* j x )(^j= 3V« (1 - VMo) I . (47) 

Obviously this can be continued, and finally we obtain the expansion for the full moment of 
inertia 

/ = jl (Mo - MoVMo + M VM VM -...)j x , (48) 
which is equivalent to the expansion of the matrix A4. in Eq. 
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